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Abstract 

The EMD algorithm, first proposed in [11], made more robust as well as more versatile in [l^, is 
a technique that aims to decompose into their building blocks functions that are the superposition of 
a (reasonably) small number of components, well separated in the time-frequency plane, each of which 
can be viewed as approximately harmonic locally, with slowly varying amplitudes and frequencies. The 
EMD has already shown its usefulness in a wide range of applications including meteorology, structural 
stability analysis, medical studies - see, e.g. [K]. On the other hand, the EMD algorithm contains 
heuristic and ad-hoc elements that make it hard to analyze mathematically. 

In this paper we describe a method that captures the flavor and philosophy of the EMD approach, 
albeit using a different approach in constructing the components. We introduce a precise mathematical 
definition for a class of functions that can be viewed as a superposition of a reasonably small number 
of approximately harmonic components, and we prove that our method does indeed succeed in decom- 
posing arbitrary functions in this class. We provide several examples, for simulated as well as real 
data. 



1 Introduction 

Time- frequency representations provide a powerful tool for the analysis of time dependent signals. They can 
give insight into the complex structure of a "multi-layered" signal consisting of several components, such 
as the different phonemes in a speech utterance, or a sonar signal and its delayed echo. There exist many 
types of time-frequency (TF) analysis algorithms; the overwhelming majority belong to either "linear" or 
"quadratic" methods. 

In "linear" methods, the signal to be analyzed is characterized by its inner products with (or correlations 
with) a pre-assigned family of templates, generated from one (or a few) basic template by simple operations. 
Examples are the windowed Fourier transform, where the family of templates is generated by translating 
and modulating a basic window function, or the wavelet transform, where the templates are obtained by 
translating and dilating the basic (or "mother") wavelet. Many linear methods, including the windowed 
Fourier transform and the wavelet transform, make it possible to reconstruct the signal from the inner 
products with templates; this reconstruction can be for the whole signal, or for parts of the signal; in the 
latter case, one typically restricts the reconstruction procedure to a subset of the TF plane. However, in all 
these methods, the family of template functions used in the method unavoidably "colors" the representation, 
and can influence the interpretation given on "reading" the TF representation in order to deduce properties 
of the signal. Moreover, the Heisenberg uncertainty principle limits the resolution that can be attained in 
the TF plane; different trade-offs can be achieved by the choice of the linear transform or the generator(s) 
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Figure 1: Examples of linear time- frequency representations. 

Top row: the signal s{t) = si{t) + S2{t) (top left) defined by si{t) = M + cos(20t) for < t < 5tt/2, (top 
middle) and S2(t) = cos (| [{t - lOf - (27r - 10)^] + 10(i - 27r)) for 2Tr < t < An (top right); Next row: 
Left:the instantaneous frequency for its two components (left) u!{t) = 20 for < (i — 10)^ < 5tt/2, and 
u}{t) = At^ + 10 for 2Tr < t < An; Middle: two examples of (the absolute value of) a continuous windowed 
Fourier transform of s{t), with a wide window (top) and a narrow window (bottom) [these are plotted 
with Matlab, with the 'jet' colormap]; Right: two examples of a continuous wavelet transform of s{t), 
with a Morlet wavelet (top) and a Haar wavelet (bottom) [plotted with 'hsv' colormap in Matlab]. The 
instantaneous frequency profile can be clearly recognized in each of these linear TF representations, but it 
is "blurred" in each case, in different ways that depend on the choice of the transform. 

for the family of templates, but none is ideal, as illustrated in Figure [l] In "quadratic" methods to build 
a TF representation, one can avoid introducing a family of templates with which the signal is "compared" 
or "measured". As a result, some features can have a crisper, "more focused" representation in the TF 
plane with quadratic methods (see Figure [2|. However, in this case, "reading" the TF representation of 
a multi-component signal is rendered more complicated by the presence of interference terms between the 
TF representations of the individual components; these interference effects also cause the "time-frequency 
density" to be negative in some parts of the TF plane. These negative parts can be removed by some 
further processing of the representation [8], at the cost of reintroducing some blur in the TF plane again. 
Reconstruction of the signal, or part of the signal, is much less straightforward for quadratic than for linear 
TF representations. 

In many practical applications, in a wide range of fields (including, e.g., medicine and engineering) one 
is faced with signals that have several components, all reasonably well localized in TF space, at different 
locations. The components are often also called "non-stationary" , in the sense that they can present jumps 
or changes in behavior, which it may be important to capture as accurately as possible. For such signals 
both the linear and quadratic methods come up short. Quadratic methods obscure the TF representation 
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Figure 2: Examples of quadratic time-frequency representations. Left: the Wigner-Ville transform 
of s{t), in which interference causes the typical Moire patterns; Middle and Right: two pseudo- Wigner- 
Ville transforms of s{t), obtained by blurring the Wigner-Ville transform slightly (middle) and somehwat 
more (right). [All three graphs plotted with 'jet' colormap in Matlab, calibrated identically] The blurring 
removes the interference patterns, at the cost of precise location in the time-frequency localization. 



with interference terms; even if these could be dealt with, reconstruction of the individual components would 
still be an additional problem. Linear methods are too rigid, or provide too blurred a picture. Figures [l] 
[2] show the artifacts that can arise in linear or quadratic TF representations when one of the components 
suddenly stops or starts. Figure |3] shows examples of components that are non harmonic, but otherwise 
perfectly reasonable as candidates for a single-component-signal, yet not well represented by standard TF 
methods, as illustrated by the lack of concentration in the time-frequency plane of the transforms of these 
signals. The Empirical Mode Decomposition (EMD) method was proposed by Norden Huang [11' as an 




Figure 3: Two examples of wave functions of the type cos[(j)(t)], with slowly varying A(t) and 4''(t), for 
which standard TF representations are not very well localized in the time-frequence plane. Left: signal. 
Middle left: windowed Fourier transform; Middle Right: Morlet wavelet transform; Right: Wigner-Ville 
function. [All TF representations plotted with 'jet' colormap in Matlab.] 

algorithm that would allow time-frequency analysis of such multicomponent signals, without the weaknesses 
sketched above, overcoming in particular artificial spectrum spread caused by sudden changes. Given a 
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signal s{t), the method decomposes it into several instrinsic mode functions (IMF): 



K 



(1.1) 



fe=i 



where each IMF is basically a function oscillating around 0, albeit not necessarily with constant frequency: 



Sk{t) = Ak{t) cosiMt)) , with Ak{t), c^'k{t) >Oyt 



(1.2) 



Essentially, each IMF is an amplitude modulated-frequency modulated (AM-FM) signal; typically, the 
change in time of Ak{t), (j)'i.{t) is much slower than the change of 4>k{t) itself, which means that locally (i.e. 
in a time interval [t — S,t + S], with S « 2'!r[(j>'f.(t)]~^) the component Sk{t) can be regarded as a harmonic 
signal with amplitude Ak{t) and frequency (j)'i.{t)- (In [TT], the conditions on an IMF are phrased as follows: 
(1) in the whole data set, the number of extrema and the number of zero crossings of Sk{t) must either 
be equal or differ at most by one; and (2) at any t, the value of a smooth envelope defined by the local 
minima of the IMF is the negative of the corresponding envelope defined by the local maxima.) After the 
decomposition of s{t) into its IMF components, the EMD algorithm proceeds to the computation of the 
"instantaneous frequency" of each component. Theoretically, this is given by Wkit) := (f'ki'^)^ practice, 
rather than a (very unstable) differentiation of the estimated ipkit), the originally proposed EMD method 
used the Hilbert transform of the Sk{t) [11]; more recently, this has been replaced by other methods [13]. 



It is obvious that every function can be written in the form (1.1) with each component as in (1.2 1. If 



s{t) is supported (or observed) in [— T, T], then the Fourier series on [—T,T] of s{t) is actually such 
a decomposition. It is also easy to see that such a decomposition is far from unique. This is simply 
illustrated by considering the following signal: 



s{t) = .25 cos([f] - 7]i) + 2.5 cos(m) + .25 cos([f] + 7]^) 



cos 



7 



) cos(m) , (1.3) 



where ^ 7, so that one can set A{t) :— 2 + cos^ [ 2 ^] ' which varies much more slowly than cos[(/)(<)] = 
cos[rit]. The interpretation of this signal is not unique: It can be regarded either as a summation of three 




Figure 4: Non- uniqueness for decomposition into IMT. This function can be considered as a single com- 
ponent, of the type A(t) cos[rit], with slowly varying amplitude, or as the sum of three components. (See 
text.) 



cosines with frequencies — 7, and 51 + 7 respectively, or as a single component with frequency Q. which 
has an amplitude A(t) that is slowly modulated. Depending on the circumstances, either interpretation 
can be the "best". In the EMD's framework, the second interpretation (single component, with slowly 
varying amplitude) is preferred when 3> 7; the EMD is typically applied when it is more "physically 
meaningful" to decompose a signal into fewer components if this can be achieved by mild variations in 
frequency and amplitude; in those circumstances, this preference is sensible. The (toy) example illustrates 
that we should not expect a universal solution to all TF decomposition problems. For certain classes of 
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functions, consisting of a (reasonably) small number of components, well separated in the TF plane, each 
of which can be viewed as approximately harmonic locally, with slowly varying amplitdes and frequencies, 
it is clear, however, that a technique that identifies these components accurately, even in the presence 
of noise, has great potential for a wide range of applications. Such a decomposition should be able to 
accommodate such mild variations within the building blocks of the decomposition. 

The EMD algorithm, first proposed in [11 , made more robust as well as more versatile in |13j (an extension 
to higher dimensions is now possible) , is such a technique. It has already shown its usefulness in a wide range 
of applications including meteorology, structural stability analysis, medical studies ~ see, e.g. [5J|SJ[TT]; a 
recent review is given in [12]. On the other hand, the EMD algorithm contains a number of heuristic and 
ad-hoc elements that make it hard to analyze mathematically its guarantees of accuracy or the limitations 
of its applicability. For instance, the EMD algorithm, uses a sifting process to construct the decomposition 



of type ( 1.1 1. In each step in this sifting process, two smooth interpolating functions are constructed (using 
cubic splines), one of the local maxima (s(i)), and one of the local minima (s(t)). From these interpolates, 
a mean curve of the signal is defined as m{t) = {s{t) + s{t))/2, which is then subtracted from the signal: 
ri{t) — s{t) —m{t). In most cases, ri is not yet a satisfactory IMF; the process is then repeated on ri again, 
etc . . .; this repeated process is called "sifting" . Sifting is done for either a fixed number of times, or until 
a certain stopping criterium is satisfied; the final remainder r„(t) is taken as the first IMF, si := r„. The 
algorithm continues with the difference between the original signal and the first IMF to extract the second 
IMF (which is the first IMF obtained from the "new starting signal" s{t) — si{t)) and so on. (Examples 
of the decomposition will be given in Section [s]) Because the sifting process relies heavly on interpolates 
of maxima and minima, the end result has some stability problems in the presence of noise, as illustrated 
in [.18 . The solution proposed in 18J addresses these issues in practice, but poses new challenges to our 
mathematical understanding. 

Attempts at a mathematical understanding of the approach and the results produced by the EMD method 
have been mostly exploratory. A systematic investigation of the performance of EMD acting on white noise 
was carried out in [T7]; it suggests that in some limit, EMD on signals that don't have structure (like 
white noise) produces a result akin to wavelet analysis. The decomposition of signals that are superpositions 
of a few cosines was studied in [16] , with interesting results. A first different type of study, more aimed 
at building a mathematical framework, is given in |14l 110] . which analyzes mathematically the limit of 
an infinite number of "sifting" operations, showing it defines a bounded operator on £oo, and studies its 
mathematical properties. 

In summary, the EMD algorithm has shown its usefulness in various applications, yet our mathematical 
understanding of it is still very sketchy. In this paper we discuss a method that captures the fiavor 
and philosophy of the EMD approach, without necessarily using the same approach in constructing the 
components. We hope this approach will provide new light in understanding of what makes EMD work, 
when it can be expected to work (and when not) and what type of precision we can expect. 



2 Synchrosqueezing Wavelet Transforms 

Synchrosqueezing was introduced in the context of analyzing auditory signals [7] ; it is a special case of real- 
location methods [TJ 12113], which aim to "sharpen" a time-frequency representation TZ(t,u!) by "allocating" 
its value to a different point (t' , uj') in the time- frequency plane, determined by the local behavior of 7?.(t, u) 
around (t, w). In the case of synchrosqueezing, one reallocates the coefficients resulting from a continuous 
wavelet transform to get a concentrated time-frequency picture, from which instantaneous frequency lines 
can be extracted. 

To motivate the idea, let us start with a purely harmonic signal. 



s{t) = Acos{ujt). 
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Take a wavelet ip that is concentrated on the positive frequency axis: ipiO = for ^ < 0. Denote by 
Ws{a, b) the continuous wavelet transform of s defined by this choice of -0. We have 



sit) a-i/2^(LJ)dt 



1 

2^ 
A 
47r 
A 
47r 



(2.1) 



a'/^{auj)e'''^. 



If 4^(0 is concentrated around ^ = loq, then Ws{a, b) will be concentrated around a ~ ujq/lu. However, the 
wavelet transform Ws{a,b) will be spread out over a region around the horizontal line a — luq/lu on the 
time-scale plane. The observation made in [7] is that although Ws{a,b) is spread out in a, its oscillatory 
behavior in b points to the original frequency lu, regardless of the value of a. 

This led to the suggestion to compute, for any (a, 6) for which Ws{a,b) ^ 0, a, candidate instantaneous 
frequency uj{a, b) by 

(2.2) 



uia,b) ^ ~iiWs{a,b)r^-^W,ia,b). 



For the purely harmonic signal s{t) = Acos{Lut), one obtains uj{a,b) ~ oj, as desired; this is illustrated 
in Figure [5] In a next step, the information from the time-scale plane is transferred to the time- frequency 






Figure 5: Left: the harmonic signal f{t) = sin(8i); Middle: the continuous wavelet transform of /; Right: 
synchrosqueezed transform of /. 

plane, according to the map (6, a) — > {b,uj{a,b)), in an operation dubbed synchrosqueezing. In [7], the 
frequency variable w and the scale variable a were "binned", i.e. Ws{a,b) was computed only at discrete 
values ttk, with ak — at-i = {Aa)k, and its synchrosqueezed transform Ts{bJ,b) was likewise determined 
only at the centers uji of the successive bins [w^ — ^AwjW^ -I- ^Aoj], with lo^ — ujg^i = Auj, by summing 
different contributions: 



T,(wf,6) = (Aw)-i VF,(afc,6)afe'/'(Aa)fc. 

ak:\^o{ai,,b)~uii\<Aui/2 



(2.3) 



The following argument shows that the signal can still be reconstructed after the synchrosqueezing. We 
have 



1 

2t: 
1 

2^ 



00 pOO 



s(e) V«) e'''^a-^dad^ 



-00 Jo 

00 pOQ 



' ^°"s(C)e'^CdC. 



(2.4) 
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Setting — 2 /q°°'0('?) we then obtain (assuming that s is real, so that s(^) = s(— 0; hence s{h) 
(47r)-i9^e [/^^ e^"? d^] ) 

r poo 

Ws{a,h) a-3/2da 



s(6) = d\t 



C-, 



In the piecewise constant approximation corresponding to the binning in a, this becomes 



s(6) « D^e 



C^i 5]t,(c.,,&) (Ac.) 



(2.5) 



(2.6) 



Remark. As defined above, (2.3 1 imphcitly assumes a Hnear scale discretization of uo. If instead logarithmic 



discretization is used, the Aw has to be made dependent on t, alternatively, one can also change the 
exponent of a from —3/2 to —1/2. 

If one chooses (as we shall do here) to continue to treat a and u> as continuous variables, without discretiza- 
tion, the analog of (2.3 1 is 



A{b) 



Ws{a,b) 5{uj{a,b)-uj)da, 



(2.7) 



where A{b) — {a; Ws{a, 6) ^ }, and uj{a, b) is as defined in (2.2) above, for (a, b) such that a e A(b). 

Remark. In practice, the determination of those (a, fe)-pairs for which Ws{a,b) = is rather unstable, 
when s has been contaminated by noise. For this reason, it is often useful to consider a threshold for 
|M^s(a, 6)1, below which uj{a, b) is not defined; this amounts to replacing A{b) by the smaller region A^{b) 
{a; \Ws{a,b)\>e}. 



One can also view synchrosqueezing as follows. For sufficiently "nice" s and ip, we have, for a S A(b), 

-id 



Lo{a,b) = --i{Ws{a,b)r' — Ws{a,b) 



/-z.-(t)a-V^V.(^)dt 
/.(i)a-i/2^(*_i^)dt 

lis + iHs)it) a-i/2 iP{*-^)dt 



(2.8) 
(2.9) 

(2.10) 

(2.11) 



where Hs denotes the Hilbert transform of s; the last equality uses that ip is supported on the positive 
frequencies only. If now s is a so-called "asymptotic signal", i.e., if 



s(t) = a{t) cos{(j){t)), 

with a'{t) ^ 1 and (l)"{t) <C 0'(i), then the Hilbert transform of s is (approximately) given by 

Hs{t) - a{t) sm{(f>{t)). 
Therefore, using the above expression of Lu{a,b), we have approximately 



, /a(.)^-(.)e-^(Oa-V^^(^)d. 

Ja{t) e**(*) a-i/2^(*-^)d< 



(2.12) 



(2.13) 



(2.14) 
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where, in the first approximation, we have omitted the term containing a'{t) as a'{t) ^ and in the 

second approximation, we have used that tp is localized around 0. 

This heuristic argument suggests that, for asymptotic signals, synchrosqueezing an appropriate wavelet 
transform will indeed give a single line on the time-frequency plane, at the value of the "instantaneous 
frequency" of a (putative) IMF. 

3 Main Result 

We define a class of functions, containing intrinsic mode type components that are well-separated, and show 
that they can be identified and characterized by means of synchrosqueezing. 

We start with the following definitions: 

Definition 3.1. [Intrinsic Mode Type Function] 

A function / : K C is said to be intrinsic-mode-type (IMT) with accuracy e > if / and A := |/| have 
the following properties; 

/(i) = A{t) e*'^(*) where A ^ {M) , (t> € C^{R) 

inicfy'it) > , 

teK 

\A'{t)\, |0"(t)| <e|0'(t)|,Vi€R 
M" := sup|(/)"(i)| < oo . 

Definition 3.2. [Superposition of Well-Separated Intrinsic Mode Components] 

A function / : M ^ C is said to be a superposition of, or to consist of, well-separated Intrinsic Mode 
Components, up to accuracy e, and with separation d if it can be written as 

m = E 

fe=i 

where all the fj. are IMT, and where moreover their respective phase functions satisfy 
<l>'kit) > <P'k-rit) , and \ct>i{t) - </.;_i(i)| > d[<t>'k{t) + <f>'k-iit)] , e R. 

Remark. It is not really necessary for the components fk to be defined on all of M. One can also suppose 

that they arc supported on intervals, supp(/fe) = supp(Afe) C [— T^.Tfe], where the different Tk need 
not be identical. In this case the various inequalities governing the definition of an IMT function or a 
superposition of well-separated IMT components must simply be restricted to the relevant intervals. For 
the inequality above on the 4>'i^{t), it may happen that some t are covered by (say) [— Tfc,Tfc] but 

not by [—Tk-i,Tk-i]; one should then replace /c — 1 by the largest £ < k for which t G [— T^,T^]; other, 
similar, changes would have to be made iftG [— Tfe_i,Tfe_i] \ [— Tfc,Tfe]. 
We omit this extra wrinkle for the sake of keeping notations manageable. 

Notation. [Class Ac,d\ 

We denote by Ae,d the set of all superpositions of well-separated IMT, up to accuracy e and with separation 
d. 

Our main result is then the following: 
Theorem 3.3. [Main result] 

Let f be a function in Ae^d, o-nd set e := e^/^. Pick a wavelet tp such that its Fourier transform 'tp 
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is supported in [1 — A, 1 + A], with A < d/(l + d), and set TZ^ = \/27r J 'ip{Q) dQ . Consider the 
continuous wavelet transform Wf{a,h) of f with respect to this wavelet, as well as the function Sf^a-ibjU) 
obtained by synchrosqueezing Wf, with threshoWe, i.e. 

Sfz{b,uj) -.^ [ Wf{a,b)S{uj-ujf{a,b))a-^/^da , 

J A,, fib) 

where A^j^b) := {a G R+ ; \Wf{a,b)\ > ?} . 

Then, provided e (and thus also 'e) is sufficiently small, the following hold: 

• \Wf{a,b)\ > e only when, for some k e {1, . . . , K} , {a,b) e Z^. := {{a,b) ; | a<p'^{b) - 1 1 < A} . 

• For each fc e {1, . . . , K}, and for each pair (a, b) e Zk for which |PV/(a, b)\ > e, we have 

\u;f{a,b) - 0^(6)1 < ?. 

• Moreover, for each k € {1, . . . , K}, there exists a constant C such that, for any & G M, 



< Ce 



This theorem basically tells us that, for / G At. a, the synchrosqueezed version Sf^i of the wavelet transform 
Wf is completely concentrated, in the (t, a;)-plane, in narrow bands around the curves uj — (j^'i^it), and 
that the restriction of S'/^e to the fc-th narrow band suffices to reconstruct, with high precision, the fc-th 
IMT component of /. Synchrosqueezing (an appropriate) wavelet transform thus provides the adaptive 
time-frequency decomposition that is the goal of Empirical Mode Decomposition. 



The proof of Theorem 3.3 relies on a number of estimates, which we demonstrate one by one, at the same 
time providing more details about what it means for e to be "sufficiently small". In the statement and 
proof of all the estimates in this section, we shall always assume that all the conditions of Theorem [3^ are 
satisfied (without repeating them), unless stated otherwise. 

The first estimate bounds the growth of the Ak, (p'f. in the neighborhood of t, in terms of the value of 
Estimate 3.4. For each fc G {1, . . . ,K}, we have 

\Ak{t+s) - Ak{t)\ < e\s\ (^|0',(i)l + and \cl,'^{t+s) - ^',{t)\ < e\s\ (^10^^)1 + ^ < kl 



Proof. 



\Akit + s) - Akit)\ = 



< 



A'^.{t + u) du 



\ A'l^it + u)\ du < e \ (j)'^{t + u)\ du 



(f>'l{t + u) du 



dt<e[ |0',(i)l kl 



The other bound is analogous. 



□ 
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The next estimate shows that, for /and -0 as given in the statement of Theorem [3^3) the wavelet transform 
Wf{a, b) is concentrated near the regions where, for some fc e {1, 2, . . . , K}, a(t>'^.{h) is close to 1. 



Estimate 3.5. 



where 



K 



K 



k=l 



fe=l 



fc=i 

within ■= J . 
Proof. We have 

K 

Wf{a,b) = Ey Mt)e'^ 



k=l 
K 



k=l V " / 



dt 



K 



It follows that 



dt 



k=l 



K 



1/2 



fe=l 



<f2 [ e\t-b\ f\cb',{b)\ + Uik\t-b\ 
k=i •' ^ 



dt 



-1/2 



k=l 



[^',{b+u)-4,',{b)]du 



t-b 



-1/2 



dt 



t^b 



dt 



K 

^^E 

/c=l 



/ \umu)\du + a^/^\M'^ j \u\^\^{u)\du 



K 



fe=i 



\\t-b\''\^',{b)\ + 



dt 



<ea'/' |/iE I'^U^')! + ^/2aE [Af^' + \A,{b)\\<P',ib)\] + ^ha^ E ^I^^Wlj 

I fe=l fc=l k=l ) 



□ 



Synchrosqueezing Wavelet Transforms - EMD 



11 



The wavelet ij: satisfies ipiC) only for 1 — A < ^ < 1 + A; it follows that 114^/(0, 6) | < e a^/^ Fi whenever 
I o.'t>'ki^) ~ 1 I > A for all fc e {1, . . . ,K}. On the other hand, we also have the following lemma: 

Lemma 3.6. For any pair {a,b) under consideration, there can be at most one k e {1, . . . ,K} for which 
\a^',ib) - 1| < A. 

Proof. Suppose that k, £ E {1,...,A'} both satisfy the condition, i.e. that |a(/)^(6) — 1| < A and 
I a4>'i{b) — 1 1 < A, with k ^ I. For the sake of definiteness, assume k> I. Since / € ^c.d, we have 

- u^) > - ^'k-iib) > d + ^'k-im > d + 0£(&)] ■ 

Combined with 

<f>'k{b) ~ 4>',{b) < a-M(l + A) - (1-A)] - 2a-' A , 
'P'kib) + m > «"M(1-A) - (1-A)] = 2a-i(l-A) , 



this gives 



A > d(l - A) 



which contradicts the condition A < d/{l + d) from Theorem 3.3 



□ 



It follows that the a, 6-plane contains K non-touching "zones", corresponding to \a(j)'i.{b) — 1| < A, 
fc € {1, . . • ,K}, separated by a "no-man's land" where |VF/(a, 6)| is small. We shall assume (see below) 
that e is sufficiently small, i.e., that for all (a, &) under consideration, 



e < a-9/4F- 



3/2 



(3.1) 



so that ea^^^Ti < e^/'^ = The upper bound in the intermediate region between the K special zones 
is then below the threshold allowed e^for the computation of LUf(a,b) used in Sf^e (see the formulation of 
Theorem 3.3 1. It follows that we will compute cjy(a, 6) only in the special zones themselves. We thus need 
to estimate dbWf{a, b) in each of these zones. 

Estimate 3.7. For k e {1,. . . , K}, and (a, 6) € M+ x M such that \ a (f>'^,{b) - 1 | < A, we have 



-idbWfiaM - V2^Afe(6)e'^'=(''Va0fe(^)^^(a'^fc(&)) 



< ea^/^Fs 



where 



K 



K 



K 



k=l 



r2 := i[Y. + ^^2«E [Mfe + \AmWkm + Afn^fc(&)i , 
fe=i fc=i 

withl'^ / \u\^\i)'{u)\du . 

Proof. The proof follows the same lines as that for Estimate |3.5| We have 



dbWf{a,b) ^ dAY. I Mt)e'*'-^'^ a'' 

^_„-3/2^ / Mt) 



i'l'lit) ^1 



dt 



t-b 



dt 



^ „ ^ f ^tl<l>i{b) + cl,',ib){t-b) + l^-''l^'t{b+ii)-^'tib)]dii] Q-3/2^/ ft b\ 

i=\ V " / 



dt 



K 



J2 IMt) - Mb)] e 



^Mt) ^-3/2 ^, 



t-b 



dt 
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By Lemma 3.6 only the term ioi £ = k survives in the sum for (a, b) such that | a(t>'^.{h) — 1 1 < A, and 
we obtain 



dbWf{a,h) - iV2TTAk{b)e"^''^'''>^(f>'k{b)^{a(j)'f^{b)) 



dbWf{a,b) - \/2^^fc(6)e''^'=(^)^^'(a0U&)) 

v a 



< J e\t-b\(^\^[ib)\ + ^M"\t-b\^ ^'(^ 



't-b 



\Mb)\ 



a 

-3/2 



dt 



t-b 



dt 



< e 



a'/'mb)\ j \u\W{u)\du + a^'^\M'; j \u\'W{u)\du 



\Mb)\' 



^|i-6n</.l.(6)| +^-\t-b\^M'^ 



<eaV2 U\cl,',{b)\ + \l'^a [M'^ + \Au{b)\\d,'^{b)\] + ^-I'^a'M'^\Ak{b)\ 



dt 



□ 



Combining Estimates 3.5 and 3.7 we find 

Estimate 3.8. Suppose that (3.1) is satisfied. For k G {\^...,K}, and {a,b) G 
\a(f>'i^{b) — 1 I < A and Wf{a, b) >'e are satisfied, we have 

\u;f{a,b) - < V^{r2 + aT,q^',{b))e'/' . 



X M such that both 



Proof. By definition, 

-idbWf{a,b) 

For convenience, let us, for this proof only, denote \/2n Ak{b) e^'^'''^^^ ^/a"^ (a 0fc(^)) by B. For the (a, 6)-pairs 
under consideration, we have then 

\ -idbWf{a,b) - (j)'i,{b)B\ < ea^^^T2 and \Wfia,b)~B\<ea^^^Ti. 

Using 7 = e^^^, it follows that 

-»9b^;(a,b) - , [B Wf{a,b)]c^',{b) 
u:Aa,b) ~ Mb) = ^^^^ + ^^^^ , 

so that 

□ 

If (see below) we impose an extra restriction on e, namely that, for all (a, 6) under consideration, and all 

fee {1,..., if}, 

e < [Fa + a </-',.(&) FJ ' , (3.2) 

then this last estimate can be simplified to 

\LOf{a,b) ~ Mb)\ (3.3) 

Next is our final estimate: 
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Estimate 3.9. Suppose that both (3.11 and (3.2 1 are satisfied, and that, in addition, for all b under 
consideration, 

e < 1/8 0^(6) + 0^(6) ]3 . (3.4) 
Let S'f^e be the synchrosqueezed wavelet transform of f , 



Sf^z{b,uj) := [ Wf{a,b) S{uj - ujf{a,b)) a-^/^ da 



Then we have, for all 6 G M, and a// /c e {1, . . . , K} 



\'^-4>i{b)\<i 



< C7. 



Proof. For later use, note first that (3.4) implies that, for all fc, £ e {1, . . . ,K}, 

d[^'^{b) + 4^',{b)\ > 21. 

We have 

Sf^a{b,uj)diLj = / / Wf{a,b)S{uj - ujf{a,b))a~^^'^dadui 

\uj-<t,'^{b)\<7 ' J\uj^4>i(b)\<7! JA.jib) 



(3.5) 



A^j(b)n{\uf(aM)-cPi{b)\<Z} 



Wf{a,b) a^^/^da . 



From Estimate 3.5 and (3.1 ) we know that |W/(a, 6) | > e only when |a(/)^(6) — 1| < A for some £ € {1, . . . , K}. 
For £ ^ fc, we have (use ( |3.5[ )) 

\cof{a,b)-cj,',ib)\ > \q^',{b) ~ <j)',{b)\ - \LUfia,b)-^',ib)\ 

> d[^^(6) + 0U&)] - e > 



which, by Estimate 3.8 implies that |a0£(6) — 1| > A. Hence 

Wf{a,b)a-^^^da 



Sf,^{b,uj)duj ^ / Wf{a,t)a-^^^da 

\uj~4,'^{b)\<Z ' JA^j{b)n{\a4>'^{b)-l\<A} 



|a0'Jh)-l|<A 



{\a4>'^{b)-l\<A}\A^jib) 



From Estimate 13.51 we then obtain 



J\uj-<t>'Jb)\<7 



< 



, , Wf{a,b)a-^/^da] - Ak{b) e^'*"'^'''^ 

'|a0'fc(fc)-l|<A / 



{|a0',(b)-l|<A}\Af.^(6) 



Wf{a,b) a~3/2 da 



< 



7^^ly2^Afc(6)e^*^W ( / 

\J\a<t>' 



|a0;(6)-l|<A) 



V^^{a(j)'kib))a-^^^da - Ak{b) e'^"^^^ 



n- 



|a0',(&)-l|<A) 



e + ea 



-3/2 



da . 
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For the first term on the right hand side, since 



J\a4>'^{b)-l\<A 



'IC-iK 

by the definition of TZ^ , and hence the first term vanishes. We thus obtain 



l'^-0'fe(b)l<? 



ML 

1- A 



1/2 



1 + A 



□ 



It is now easy to see that all the Estimates together provide a complete proof for Theorem |3.3[ 

Remark. We have three different conditions on e, namely (3.1), (3.2 1 and (3.4 1. The auxiliary quantities 



in these inequalities depend on a and b, and the conditions should be satisfied for all (a, 6)-pairs under 
consideration. This is not really a problem: the dependence on b is via the quantities Ak{b) and 4'kib)j 
and it is reasonable to assume these are uniformly bounded above and below (away from zero); because 
the bounds on the 0jj(&) translate into bounds on a, we can likewise safely assume that a is bounded 
above as well as below (away from zero). Note that the different terms can be traded off in many other 
ways than what is done here; no effort has been made to optimize the bounds, and they can surely be 
improved. The focus here was not on optimizing the constants, but on proving that, if the rate of change 
(in time) of the Ak{b) and the (t)'f.{b) is small, compared with the rate of change of the (j>k{b) themselves, 
then synchrosqueezing will identify both the "instantaneous frequencies" and their amplitudes. 
In the statement of Theorem |3.3| we required the wavelet ip to have a compactly supported Fourier 
transform. This is not absolutely necessary; it was done here for convenience in the proof. If is not 
compactly supported, then extra terms occur in many of the estimates, taking into account the decay of 
i^iC) as ^ ^ oo or C ^ 0; these can be handled in ways similar to what we saw above, at the cost of 
significantly lengthening the computations without making a conceptual difference. 



4 A Variational Approach 

The construction and estimates in the previous section can also be interpreted in a variational framework. 

Let us go back to the notion of "instantaneous frequency." Consider a signal s{t) that is a sum of IMT 
components Si(t): 

N N 

s{t) = s^(^) = ^'(^) cos(0.(t)), (4.1) 

i=l i=l 

with the additional constraints that (p'^ (t) and (j)'^ {t) for i j are "well separated" , so that it is reasonable 
to consider the Si as individual components. According to the philosophy of EMD, the instantaneous 
frequency at time t, for the i-th component, is then given by uJi(t) — 4>[{t). 

How could we use this to build a time-frequency representation for si If we restrict ourselves to a small 
window in time around T, of the type [T - M,T + M], with At w 2'k/4>[{T), then (by its IMT nature) the 
i-th component can be written (approximately) as 

S^ii)\[T-At^T+At] ~ MT) COS [0,(r) -|- cP'^{T) {t - T)] , 
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which is essentially a truncated Taylor expansion in which terms of size 0{A'^{T)), 0{(l)['{T)) have been 
neglected. Introducing u!i{T) = 4i[{T), and the phase fi{T) := ^^(r) — uji{T)T we can rewrite this as 



[T-At,T+At] 



A,{T) COS [oj^{T)t + ^,{T)]. 



This signal has a time- frequency representation, as a bivariate "function" of time and frequency, given by 
(for t e [T- Ai,T + At]) 



F,{t,u) = A,{T) cosM + (^,(r)](5(w-w,(T)), 



(4.2) 



where 5 is the Dirac-delta measure. The time-frequency representation for the full signal s = X^i^i ^^^^ 
in the neighborhood oi t = T, would then be 



N 



F{t,uj)^J2MT) cos[ujt + ^,{T)]6{u;-cu,{T)). 



(4.3) 



Integrating over ui, in the neighborhood oi t — T, leads to s{t) « / F{t,uj) duj. 



All this becomes even simpler if we introduce the "complex form" of the time-frequency representation: 
for t near T, we have F{t,u) = 'E^^i^ji^) exp{iLut) S{uj - 0j(r)), with Aj{T) ^ Aj{T) exp[i^j(T)]; 
integration over ui now leads to 



F{t,iu)du 



N 



exp[iu;,{T)t + i^,{T)] 



sit). 



(4.4) 



Note that, because of the presence of the ^- measure, and under the assumption that the components remain 
separated, the time-frequency function F(t,uj) satisfies the equation 



dtF{t,uj) = iLuF{t,uj). 



(4.5) 



To get a representation over a longer time-interval, the small pieces described above have to be knitted 
together. One way of doing this is to set F{t,uj) — J2f=i exp{iujt) 5{uj — ujj{t)). This more globally 

defined F(t,Lu) is still supported on the N curves given by = Wj{t), corresponding to the instantaneous 
frequency "profile" of the different components. The complex "amplitudes" Aj(t) are given by Aj(t) = 
Aj(t) exp[iipj{t)], where, to determine the phases exp[iipj(t)], it suffices to know them at one time to- We 
have indeed 

= ^ m) ^Am = m ^At) ^'m = -^(*)*; 

since the uij (t) are known (they are encoded in the support of F) , we can compute the ipj (t) by using 



ifi^it) 



to 



w^(r)rdT -I- (pj{to)- 



Moreover, (4.5 1 still holds (in the sense of distributions) up to terms of size 0(^-(r)), 0((/)-'(r)), since 

N 

dtF{t,u) = { [A^t) - iujAt)tAj{t) + iLuA^it)] e"^* 5{uj - uj^it)) + Aj{t) e"^* uj^t) - t^j(t))} 



N 



= iu:Y,A,{t)e^-U{u-uj,{t)) + 0{^,{T), c^'l{T)) 
= iL0F{t.^) + O{A[{T),cj>'l{T)). 
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This suggests modeling the adaptive time-frequency decomposition as a variational problem in which one 
seeks to minimize 



me 



F{t,Lu)dLO 



-sit) 



dt + fi // \dtF{t,uj)-iujF{t,Lu)\ dtduj 



(4.6) 



to which extra terms could be added, such as, 7 // \F{t,uj)\'^ dtdcu (corresponding to the constraint that 
F e _L^(]R^)), or A J [[/ |F(t,ti;)| dw]^ dt (corresponding to a sparsity constraint in lo for each value of t). 
Using estimates similar to those in Section [3] one can prove that if s G Ae,d, then its synchrosqueezed 
wavelet transform Ss,7(b,Lo) is close to the minimizer of (4.6). Because the estimates and techniques of 
proof are essentially the same as in Section [3] we don't give the details of this analysis here. 

Note that wavelets or wavelet transforms play no role in the variational functional - this fits with our 
numerical observation that although the wavelet transform itself of s is definitely influenced by the choice 
of ip, the dependence on ^ is (almost) completely removed when one considers the synchrosqueezed wavelet 
transform, at least for signals in A^^d- 



5 Numerical Results 

In this section we illustrate the effectiveness of synchrosqueezed wavelet transforms on several examples. 
For all the examples in this Section, synchrosqueezing was carried out starting from a Morlet Wavelet 
transform; other wavelets that axe well localized in frequency give similar results. 



5.1 Instantaneous Frequency Profiles for Synthesized data 

We start by revisiting the toy signal of Figures [l] and [2] in the Introduction. Figure [6] shows the result 
of synchrosqueezing the wavelet transform of this toy signal. We next explore the tolerance to noise of 






iime 



time 



time 



Figure 6: Revisiting the toy example from the Introduction. Left: the toy signal used for Figures 
[1] and [2] Middle: its instantaneous frequency; Right: the result of synchrosqueezing for this signal. The 
"extra" component at very low frequency is due to the signal's not being centered around 0. 

synchrosqueezed wavelet transforms. We denote by X{t) a white noise with zero mean and variance = 1. 
The Signal-to-Noise Ratio (SNR) (measured in dB), will be defined (as usual) by 



SNR [dB] ^ 10 log 



10 



Var/ 



where / is the noiseless signal. Figure [7] shows the results of applying our algorithm to a signal consisting 
of one single chirp function /(t) = cos(8i + t^), without noise (i.e. the signal is just /), with some noise 
(the signal is / + X, SNR= -3.00 dB), and with more noise(/ + 3X, SNR= -12.55dB). Despite the high 
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Figure 7: Top row: left: single chirp signal without noise; middle: its continuous wavelet transform, and 
right: the synchrosqueezed transform. Middle row: same, after white noise with SNR of —3.00 dB was 
added to the chirp signal. Lower row: same, now with white noise with SNR of —12.55 dB. 




5 10 
time 



Figure 8: Far left: the instantaneous frequencies, cj(i) — 2t + \ — sin(i) and 8, of the the two IMT 
components of the crossover signal f{t) = cos[t^ + t + cos(t)] + cos(8f); Middle left: plot of f{t) with 
no noise added; Middle: Synchrosqueezed wavelet transforms of noiseless f{t); Middle right: /(<)+noise 
(corresponding to SNR of 6.45 dB); Far right: synchrosqueezed wavelet transforms of /(t)+noise. 



noise levels, the synchrosqueezing algorithm identifies the component with reasonable accuracy. Figure 10 
below shows the instantaneous frequency curve extracted from these synchrosqueezed transforms, for the 
three cases. 

Finally we try out a "crossover signal", that is, a signal composed of two components with instantaneous 
frequency trajectories that intersect; in our example f{t) = cos(i^ + i + cos(i)) + cos(8t). Figure [s] shows 
the signals f(i) and /(t) +0.5X(i), together with their synchrosqueezed wavelet transforms, as well as the 
"ideal" frequency profile given by the instantaneous frequencies of the two components of /. 
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5.2 Extracting Individual Components from Synthesized data 





5 10 
time 





5 10 
time 




Figure 9: Comparing the decomposition into components si{t) and S2{t) of the crossover signal 

s{t) — si{t) + S2{t) — cos(8<) + cos[t^ + i+cos(i)] Rows: Noise-free situation in the first two rows; in the last 
two rows noise with SNR — 6.45 dB was added to the mixed signal. In each case, Si is in the top row, and 
S2 underneath. Columns: Far left: true Sj{t) j — 1,2; Middle left: Zone marked on the synchrosqueezed 
transform for reconstruction of the component; Center: part of the synchrosqueezed transform singled out 
for the reconstruction of a putative sj; Middle Right: the corresponding candidate Sj(t) according to the 
synchrosqueezed transform (plotted in blue over the original Sj, in red); Far right: candidate Sj according 
to EMD in the noiseless case, EEMD in the noisy case. 



In many applications listed in [SJinillllllllllH], the desired end goal is the instantaneous frequency trajectory 
or profile for the different components. When this is the case, the result of the synchrosqueezed wavelet 
transform, as illustrated in the preceding subsection, provides a solution. 

In other applications, however, one may wish to consider the individual components themselves. These are 
obtained as an intermediary result, before extracting instantaneous frequency profiles, in the EMD and 
EEMD approaches. With synchrosqueezing, they can be obtained in an additional step after the frequency 
profiles have been determined. 

Recall that, like most linear TF representations, the wavelet transform comes equipped with reconstruction 
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formulas, 



fit) = a 



OO pOO 



V, / / Wf{a,b)a-^/'^^ i^- — ^) dadb 

-OO Jo 



as well as /(t) = C'^ f Wf{a,t)a~^^^da^C'^me [ 
Jo . Jo 



Tf(t, uj) duj 



(5.1) 
(5.2) 



where C^, are constants depending only on ?/;. For the signals of interest to us here, the synchrosqueezed 
representation has, as illustrated in the preceding subsection, well localized zones of concentration. One 
can use these to select the zone corresponding to one component, and then integrate, in the reconstruction 
formulas, over the corresponding subregion of the integration domain. In practice, it turns out that the 



best results are obtained by using the reconstruction formula (5.1 ) 



We illustrate this with the crossover example from the previous subsection: Figure [9] shows the zones 
selected in the synchrosqueezed transform plane as well as the corresponding reconstructed components. 
This figure also shows the components obtained for these signals by EMD for the clean case, and by EEMD 
(more robust to noise than EMD) for the noisy case; for this type of signal, the synchrosqueezed transform 
proposed here seems to give a more easily interpretable result. 

Once the individual components are extracted, one can use them to get a numerical estimate for the 
variation in time of the instantaneous frequencies of the different components. To illustrate this, we revisit 
the chirp signal from the previous subsection. Figure [lO] shows the frequency curves obtained by the 
synchrosqueezing approach; they are fairly robust with respect to noise. 






Figure 10: Instantaneous frequency curves extracted from the synchrosqueezed representa- 
tion Left: the instantaneous frequency of the clean single chirp signal estimated by synchrosqueezing; 
Middle: the instantaneous frequency of the fairly noisy single chirp signal (SNR=— S.OOdB) estimated by 
synchrosqueezing. Right:the instantaneous frequency of the noisy single chirp signal (SNR=— 12.55dB) 
estimated by synchrosqueezing. 



5.3 Applying the synchrosqueezed transform to some real data 

So far, all the examples shown concerned toy models or synthesized data. In the subsection we illustrate 
the result on some real data sets, of medical origin. 

5. 3. A Surface Electromyography Data 

In this application, we "clean up" surface electromyography (sEMG) data acquired from a healthy young 
female with a portable system (QuickAmp). The sEMG electrodes, with sensors of high-purity sintered 
Ag/AgCl, were placed on the triceps. The signal was measured for 608 seconds, sampled at 500Hz. During 
the data acquisition, the subject flexed/extended her elbow, according to a protocol in which instructions 
to flex the right or left elbow were given to the subject at not completely regular time intervals; the subject 
did not know prior to each instruction which elbow she would be told to flex, and the sequence of left/right 
choices was random. The raw sEMG data si{t) and Sr{t) are shown in the left column in Figure 11 
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Figure 11: Left column: The surface electromyography data described in the text. The erratic (and 
different) drifts present in these signals are pretty typical. The very short-lived peaks correspond to the 
elbow-flexes by the subject. Note that the peaks have, on average, the same amplitudes for the left and 
right triceps; they look larger for the right triceps only because of a difference in scale. Middle column: 
The synchrosqueezing transforms of the surface electromyography signals. Coarse scales are near the top 
of these diagrams, finer scales are shown lower. The peaks are clearly marked over a range of fine scales. 
Right column: The red curve give the signals "sans drift", reconstructed by deleting the low frequency 
region in the synchrosqueezed transforms; comparison with the original sEMG signals (in blue) shows the 
R peaks are as sharp as in the original signals, and at precisely the same location. 



The middle column in Figure 11 shows the results Ws£{t,uj), Wsr{t,uj) of our synchrosqueezing algorithm 
applied to the two sEMG data sets. We used an implementation in which each dyadic scale interval 
(a e [2'^,2'^+^]) was divided into 32 equi-log-spaced bins. 

The original surface electromyography signals show an erratic drift, which medical researchers wish to 
remove without losing any sharpness in the peaks. To achieve this, we identified the low frequency com- 
ponents (i.e. the dominant components at frequencies below 1 Hz) in the signal, and removed them 
before reconstruction. More, precisely, we defined Si(i) = X]4>4i „„t ?)i with frequency cut-off 

cut-off (^) = ^i{t) + ^Oi where uJi{t) was the dominant component for signal i {i = £ or r) near 1 Hz with 
the highest frequency, and ujq a small constant offset. The right column in Figure [Tl] shows the results; 
the erratic drift has disappeared and the peaks are well preserved. Note that a similar result can also be 
obtained straightforwardly from a wavelet transform without any synchrosqueezing. 

5.3.B. Electrocardiogram Data 

In this application we use synchrosqueezing to extract the heart rate variability (HRV) from a real elec- 
trocardiogram (ECG) signal. The data was acquired from a resting healthy male with a portable EGG 
machine at sampling rate lOOOHz for 600 seconds. The samples were quantized at 12 bits across ±10 mV. 



The raw EGG data e{t) are (partially) shown in the left half of Figure 12 the right half of the figure shows 
a blow-up of 20 seconds of the same signal. 



The strong, fairly regular spikes in the EGG are called the R peaks; the heart rate variability (HRV) time 
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series is defined as the sequence of time differences between consecutive R peaks. (The interval between 
two consecutive R peaks is also called the RR interval.) The HRV is important for both clinical and basic 
research; it reflects the physiological dynamics and state of health of the subject. (See, e.g., [15J for clinical 
guidelines pertaining to the HRV, and [2] recent advances made in research.) The HRV can be viewed as 
a succession of snapshots of an averaged version of the instantaneous heart rate. 



The left half of Figure 13 shows the synchrosqueezed transform Te{uj,t) of e{t); in this case we used an 
implementation in which each dyadic scale interval (a £ [2'^,2'^+^]) was divided into 128 equi-log-spaced 
bins. The synchrosqueezed transform Te{LU,t) has a dominant line c(t) near 1.2Hz, the support of which 



can be parametrized as {(t,uJc{t)); t £ [0, 80sec] }. The right half of Figure 13 tracks the dependence on 
t of u!c{t). This figure also plot a (piecewise constant) function f{t) that tracks the HRV time series and 
that is computed as follows: if t lies between t — i and i^+i, the locations in time for the i-th and {i + l)-st 
R peaks, then f{t) = [ti+i — ti]^^ . The plot for Lo{t) and /(t) are clearly highly correlated. 




300 320 3-40 3SO 3SO AGO 300 320 3AO 3SO 330 400 

time time 



Figure 13: Left: The synchrosqueezing transforms of the electrocardiogram signals given in Fig. [12] 
Right: The blue curve shows the "instantaneous heart rate" a;(t) computed by tracking the support of the 
dominant curve in the synchrosqueezed transform T^] the red curve is the (piecewise constant) inverse of 
the successive RR^ 
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